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Abstract 

Positron  Survival  In  Type  II  Supernovae 
by 

Steven  J.  Sturner 

In  this  work  I  investigate  the  possibility  of  Type  II  supemovae  being  the  origin  for 
positrons  producing  observed  annihilation  radiation  observed  toward  the  Galactic  center. 
It  was  my  contention  that  the  decay  of  56Co  coupled  with  falling  densities  would  allow 
for  the  production  and  extended  existence  of  positrons  in  the  supernova  outflow. 
Supernova  1987A  has  prompted  many  people  to  construct  models  of  supernova  outflow. 
I  use  the  results  of  two  existing  models  as  the  initial  conditions  in  my  models.  I  have 
created  both  an  analytic  and  a  computer  model  for  the  survival  of  positrons.  These 
models  show  that  while  Type  II  supemovae  fall  short  of  the  needed  production  of 
surviving  positrons,  the  lower  densities  existing  in  Type  I  supemovae  may  be  a  more 
promising  source. 
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I.  Introduction 

/ 

Since  the  discovery  of  Supernova  1987A  (a  Type  II  supernova)  in  February  of  1987 
there  has  been  a  flurry  of  activity  in  the  field  of  supernova  modelling.  The  fact  that  the 
supernova  occurred  in  the  nearby  Large  Magellanic  Cloud  and  that  the  progenitor  star  has 
almost  universally  been  identified,  as  Sk  -69  202  (Woosley  1988;  Woosley  (in  press); 
Nomoto  et  al.  1988;  Shigeyama,  'Nomoto,  and  Hashimoto  1988)  has  allowed  great 
progress  to  be  made  in  this  area.  Pnor  to  exploding,  Sk  -69  202  was  seen  as  a  B3-I  blue 
supergiant  with  a  surface  temperature  of  approximately  16,000  K  and  a  bolometric 
magnitude  of  -7.5  to  -8.2.  A  blue  supergiant  becoming  a  Type  II  supernova  shocked 
many  people.  It  has  been  suggested  that  Sk  -69  202  was  once  a  red  supergiant  but  it 

evolved  back  to  the  blue.  Models  have  used  low  stellar  metallicity  and  mass  loss  to  get 

A  ?  -  -  -  i  .  ■  ,  /  ' 

such  a  star  to  evolve  back  before  exploding  (Woosley  1988;  Shigeyama,  Nomoto,  and 
Hashimoto  1988;  Nomoto  et  al.  1988).  As  a  result  of  the  proximity  of  SN  1987 A,  there 
is  a  wealth  of  data  on  the  light  curve.  Details  of  the  middle  light  curve  have  put  limits  on 
the  amount  of  radioactive  56Ni  (and  thus  the  56Co)  that  was  produced  in  the  explosive 
silicon  burning  that  took  place  as  the  shock  propagated  outward.  Details  of  the  middle 
light  curve  together  with  the  early  emergence  of  x-  and  y-rays  indicate  that  this  nickel  was 
mixed  radially  outward  as  well  (Woosley  1988;  Woosley  (in  press);  Pinto  and  Woosley 
1988;  Shigeyama,  Nomoto,  and  Hashimoto  1988;  Nomoto  et  al.  1988).  This  radial 
mixing  was  driven  by  Rayleigh-Taylor  instabilities  caused  by  the  deposition  of  y-ray 
energy  in  the  region  of  nickel  production. 

These  new  models  allow  many  calculations  to  be  made  on  the  effects  of  supemovae 
on  their  galactic  neighborhood.  One  of  these  will  be  investigated  in  this  thesis;  the 
possible  contribution  of  positrons  from  a  Type  II  supernova  to  the  observed  Galactic  5 1 1 
keV  y-ray  line.  Observations  of  celestial  annihilation  radiation  was  pioneered  by  Haymes 
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et  al.  at  Rice  University  in  the  late  60's  and  early  70's.  They  found  a  source  of  possible 
annihilation  radiation  towards  the  Galactic  center.  However,  the  source  of  the  observed 
line  was  never  pinpointed  as  positron  annihilation  because  of  inadequate  energy 
resolution.  The  measured  flux  was  1.8xl0‘3  ycnr2  s4  (if  inteipretted  as  a  point  source) 
when  using  a  detector  with  an  aperture  of  24*  FWHM  in  1970  and  1971.  Later,  in  1974, 
a  detector  with  an  aperture  of  15*  FWHM  was  used  and  a  flux  of  0.80xl0*3  ycnr2  s*1 
was  observed.  Thus,  evidence  for  a  variable  or  diffuse  source  was  seen  from  the 
beginning.  Several  groups  have  continued  this  work.  In  1977  Bell/Sandia  Labs  were 
able  to  determine  that  the  line  was  due  to  annihilation  of  positrons  because  of  their 
improved  energy  resolution.  They  observed  a  flux  of  1.22xl0*3y  cm*2  s4  towards  the 
Galactic  center  with  a  15*  FWHM  aperture.  Also  in  1977  the  Centre  d'Emde  Spatiale  des 
Rayonnements  measured  a  flux  of  4.2  or  3.2  xlO*3  y  cm-2  s4  (depending  on  reduction 
method)  using  an  aperture  of  50*  FWHM.  Comparing  their  observations  with  earlier 
ones,  they  concluded  that  the  observed  flux  depended  on  the  instrument  aperture.  This 
high  flux  was  confirmed  by  a  group  at  the  University  of  New  Hampshire  with  a  detector 
with  an  aperture  of  100*  FWHM  in  1977.  Riegler  et  al.  made  observations  using  the 
HEAO-3  satellite  (35*  FWHM)  in  1979  and  1980.  They  found  fluxes  of  1.85  and  0.65 
xl 0'3  ycnr2  s4  respectively  (Share  et  al.  1988).  Ramaty  and  Lingenfelter  saw  this  as 
evidence  of  a  variable  compact  source  that  turned  on  in  1977  and  off  in  1979  plus  a 
diffuse  source  (Ramaty  and  Lingenfelter  1987).  Share  et  al.  made  observations  using  the 
130*  FWHM  y-ray  spectrometer  on  the  Solar  Maximum  Mission  satellite  from  October  to 
February  yearly  from  1980-81  to  1985-86.  They  report  a  flux  of  about  2.1xl0-3  y  cm*2 
s4  with  less  than  a  30%  yearly  variation.  This  was  significantly  larger  than  observations 
made  with  smaller  apertures  within  that  same  time  period  (see  figure  1).  Thus  Share  et 
al.  conclude  that  the  best  explanation  is  a  diffuse  source  of  annihilation  radiation  and  the 
case  for  a  variable  source  was  weakened  because  of  the  short  period  for  variability 
needed  (Share  et  al.  1988). 


History  of  Galactic  center  511  keV  measurements.  Note  the  general  trend  of  larger 
apertures  obtaining  larger  fluxes  and  the  suggestion  of  a  variable  source  turning  on  in 
1977  and  off  in  1979  (Share  et  al.  1988). 
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Clayton  (1973)  suggested  that  the  source  of  these  positrons  should  be  either  56Co 
(daughter  of  56Ni)  or  ^Ti.  He  predicted  that  both  of  these  nuclei  should  be  produced 
during  explosive  nucleosynthesis.  The  production  of  this  nickel  has  been  proven  by  the 
light  curve  of  SN  1987 A.  He  concluded  that  both  were  potentially  adequate  to  account 
for  the  observed  feature. 

In  this  thesis  I  discuss  the  escape  of  positrons  produced  by  the  decay  of  cobalt  to 
iron  in  a  Type  II  supernova.  I  decided  not  to  look  at  whether  positrons  could  be  "shot" 
out  of  a  supernova,  but  whether  a  positron  could  survive  for  an  extended  period  of  time 
within  the  expanding  gas.  My  motivations  for  this  were  two-fold: 

1)  The  uncertain  structure  of  the  magnetic  fields  within  the  supernova  made  it  uncertain 
that  a  charged  particle  had  any  chance  to  escape. 

2)  The  velocity  profile  of  an  expanding  supernova  provided  an  environment  where  the 
density  of  electrons  would  decrease  rapidly  (section  IV). 

The  uncertain  nature  of  the  magnetic  fields  made  them  hard  to  incorporate  into  my  model. 
The  condition  that  the  positrons  were  required  to  remain  within  the  shell  in  which  they 
were  produced  implied  the  existence  of  these  fields. 

The  organization  of  this  thesis  is  as  follows.  In  sections  II- V,  I  discuss  the 
processes  involved  in  the  production,  thennalization  and  annihilation  of  positrons.  In 
section  VI,  I  discuss  published  models  of  SN  1987A  from  which  I  determined  the  initial 
conditions  for  my  models.  Sections  VII  and  VIII  contain  these  models  of  positron 
survival.  There  is  both  an  analytic  and  a  computer  model 
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n.  The  8*  Decay  of  56Co  to  56Fe 

The  decay  of  cobalt  to  iron  is  an  important  process  because  the  energy  liberated 
powers  a  supernova  during  the  middle  part  of  its  light  curve.  We  are  particularly 
interested  in  the  emission  of  positrons  from  this  reacdon  as  a  possible  source  for  the 
diffuse  component  of  the  background  511  keV  y-rays.  The  cobalt  is  formed  by  the  decay 
of  56Ni.  This  isotope  of  nickel  has  a  half-life  of  6. 1  days;  therefore,  it  all  has  essentially 
turned  to  cobalt  at  times  of  interest  to  us.  I  have  therefore  assumed  that  56Co  was  the 
explosive  nucleosynthesis  product.  The  half-life  of  56Co  is  78.8  days  (see  figure  2).  A 
positron  is  emitted  when  the  cobalt  decays  to  an  excited  state  of  iron.  This  route  is  taken 
20%  of  the  time. 


Co  56 


______  Fe  56 

The  Q-value  for  electron  capture  is  4.57  Mev  (Lederer  et  al.  1967).  From  this,  one 
can  determine  that  the  maximum  kinetic  energy  of  an  emitted  positron  is  1.46  MeV.  The 
positron  energy  distribution  function  is  of  the  form  (Whaling  1960) : 

N  =  kpw(w0  -  w)2 


where  k  =  a  normalization  constant. 
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p  =  momentum  of  the  positron, 

w  =  total  energy  of  the  positron  =  KE  +  me2, 

w0=  maximum  total  energy  of  the  positron  =  KEmax  +  me2. 

Thus  N  =  kpw[  KEmM  -  KE  ]2  .  This  can  be  simplified  so  the  only  variable  is  the  kinetic 
energy: 

1  T 

l.  Tr  ,i2r  „_i2 

N  =  |KE>KE  +  2mc2  ■  [KE  +  me2  j  [KE,,,,,  -  KE  i  .  (0-!) 


To  numerically  follow  the  thermalization  of  these  positrons  I  will  find  it  convenient  to 
approximate  the  spectrum  as  a  weighted  sequence  of  monoenergetic  emissions.  Breaking 
the  energy  range  (0-1.46  MeV)  into  ten  intervals,  I  calculate  the  fraction  of  positrons 
emitted  in  each  interval  by  finding  the  area  under  a  plot  of  equation  (II- 1)  (see  figure  3)  to 
be: 


Table  II-I:  Fraction  of  Positrons  in  Kinetic  Energy  Interval 


Interval  (MeV) 

Midpoint  Energy  of  Interval  (MeV) 

Fraction  of  Total  Positrons 

0.000  -  0.146 

0.073 

0.061 

0.146  -  0.292 

0.219 

0.117 

0.293  -  0.438 

0.365 

0.152 

0.439  -  0.584 

0.511 

0.165 

0.585  -  0.730 

0.657 

0.158 

0.731  -  0.876 

0.803 

0.142 

0.877  -  1.022 

0.949 

0.104 

1.023  -  1.168 

1.095 

0.066 

1.169  -  1.314 

1.241 

0.030 

1.315  -  1.460 

1.387 

0.005 

This  fraction  will  be  used  to  weight  quantities  when  in  determining  a  mean  value  with 
respect  to  energy. 


9 


ITT.  Positron  Energy  Loss  in  a  Partially  Ionized  Medium 

Once  a  positron  is  emitted,  it  loses  energy  in  a  variety  of  ways  as  it  passes  through 
matter.  It  will  lose  energy  to  the  un-ionized  component  by  inelastic  collisions  with  the 
atomic  electrons.  In  this  process  the  electron  will  be  removed  from  the  atom  or  it  will 
move  to  an  excited  state.  The  positrons  will  lose  energy  to  the  ionized  component  by 
interacting  with  the  free  electrons.  The  positrons  will  also  radiate  energy  via 
Bremsstrahlung  when  interacting  with  nuclei.  I  derive  expressions  for  these  in  order  to 
calculate  thermalization  times  for  positrons  emitted  by  56Co. 

nia)  Coulomb  collisions  with  bound  electrons 

The  energy  loss  rate  for  inelastic  coulomb  collisions  with  bound  electrons  has  been 
calculated  by  many  people.  The  general  form  of  this  loss  rate  can  be  gotten  by  a  simple 
classical  calculation  (Whaling  1960).  Consider  an  ion  of  mass  M,  charge  ze,  and 
velocity  v.  Let  b  be  the  impact  parameter. 

M,  ze 


The  momentum  gained  by  the  electron  is: 


f  F  dt  where  F  =  -  ze-y 
J  x2  +  b2 
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If  the  electron  is  assumed  to  move  very  little  during  the  encounter  this  can  be  simplified 
to: 


where  F  =  F  sin©  = 


With  the  substitution  dx  =  v  dt. 


F,  dx 


If  the  electron  is  regarded  as  free,  then  it  acquires  kinetic  energy: 


AE  = 


(Apf 

2m 


The  electron  may  be  quasi-free  if  the  Coulomb  force  from  the  passing  ion  greatly  exceeds 
that  holding  the  electron  in  its  orbit.  Therefore  the  amount  of  energy  transferred  from  the 
ion  to  a  bound  electron  is: 


AE  = 


2z2e4 

b2v2m 


If  there  are  NZ  electrons  per  cubic  centimeter,  then  the  ion  will  encounter  2icb  db  dx  NZ 
electrons  with  an  impact  parameter  between  b  and  b+db  in  traveling  a  distance  dx.  Thus 
the  energy  lost,  dE,  will  be  27tb  db  dx  NZ  AE.  The  total  amount  of  energy  lost  in 
traveling  a  distance  dx  is  found  by  integrating  over  all  impact  parameters: 


dx 


mv2 


hljtl 

f  \ 

U 

f  db  _  4:tz2e4NZ 

^mtx 

J  b  "  mv2  1 

1 

0  .  * 
min 

v  y 
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if  z=l  (as  for  positrons) 


dE  _ 
dx 


(  \ 

(  \ 

(  \ 

4jte4  I 

NZ 

^max 

lmc2  J1 

p2 

v  ; 

In 

b  . 

1 mm ) 

(HI-1) 


The  values  of  and  are  estimated  in  different  ways.  One  example  is: 


b  . 

min 


=  max 


Ze2  h 
j^ymv2  ’  V™ 


1 

I 

J 


.  •  ^ 
Dmax  —  ^  . 


where  y  is  the  relativistic  correction  factor  and  o>  is  a  characteristic  atomic  frequency  of 
motion  (Jackson  1975).  Physically,  here  the  minimum  impact  parameter  is  either  the 
impact  parameter  where  the  maximum  allowable  energy  (2y2mv2)  is  transferred  to  the 
bound  electron  using  the  above  form  for  AE  or  the  minimum  extent  of  the  electron's 
wave  packet  The  maximum  impact  parameter  is  gotten  by  setting  the  collision  time  equal 
to  the  orbital  period.  If  the  collision  time  is  greater  than  the  orbital  period  there  is  no  net 
transfer  of  energy  to  the  atomic  electron. 

More  in  depth  calculations  yield  various  results.  Bethe  derived  a  quantum 
mechanical  formula  for  the  energy  loss  of  an  electron  owing  to  ionization  of  bound 
electrons  as  it  passes  through  matter  with  relativistic  corrections  and  corrections  due  to 
the  electron  spin  (Heitler  1934;  Wu  1960): 

dE  _  27Ce4NZbound 
dx  mv2 


In 

y^mv^E  ^ 
-2 

+  A 

l  21  J 

-r 


pel  ’ 


(m-2) 


where  I  is  the  mean  ionization  potential  (»11.5  eV  x  Z), 
NZbound  is  the  density  of  atomic  electrons, 

P  =  v/c, 


-1 
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*-k  1-  cl 


A  =  ( 1  -  P2  ) - (iTTp1  - 1  +  P2  )n2  +  i(  1 .  yr? ) , 


and  r  j  is  a  correction  term  resulting  from  dielectric  properties  of  the  medium. 


Equation  (IH-2)  can  be  transformed  into: 


dE  _  47te4 


NZbo, 


r  r 


^  ini  pi  E±mci  LeJ  |-iin2+lA  }-r 

P  '{  L  l  1  Mc  >  J  J 


(m-3) 


This  gives  a  factor  multiplying  the  logarithm  that  is  the  same  as  in  the  classical  equation 
(HI-1).  In  order  for  equation  (III-2)  to  apply  to  positrons  I  believe  that  the  argument  of 
the  logarithm  should  be  increased  by  a  factor  of  two  in  equation  (DI-2).  For  the  case  of 
an  incoming  electron  arguments  were  made  that  the  maximum  energy  transfer  would  be 
1/4  mv2  instead  of  1/2  mv2  because  it  is  impossible  to  distinguish  whether  the  detected 
outgoing  electron  is  the  target  or  the  projectile.  But  a  positron  and  electron  are 
distinguishable.  This  would  eliminate  the  -1/2  In  2  in  equation  (DI-3).  The  details  of  the 
calculation  were  not  provided  making  it  impossible  to  tell  for  certain.  Since  it  is  not  the 
purpose  of  this  thesis  to  investigate  details  of  energy  loss,  I  will  assume  equation  (HI-3) 
is  sufficiently  correct  for  positrons. 

Another  version  of  this  formula  was  listed  for  positrons  or  electrons  by  M.  Zombeck 
in  Special  Report  386  of  the  Smithsonian  Astrophysical  Observatory  (Sect  13,  p.  23, 
1980): 

(  F  .  ,  .1  1 

=  ^  In  [  E+mC2  V  _E_  |  I  .  ifl2  ^ 

dx  «s  p2  s—  j  LPl— |mc2  J  j  Jr 
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f 


4  tre4 


me 


NZ.. 


ixjund 


V 


r  r 

•I  In 

l  L 


■  1  1 


p 


E+mc2 


me 


-ip2  K 

J  J 


(in-4) 


A  “ 

where  r0  =  .  =  the  classical  electron  radius  , 

me2 

and  as  before,  I=the  mean  ionization  potential,  (3=v/c,  and  NZbo^sthe  density 
of  atomic  electrons. 


A  third  version  of  this  formula  is  given  for  electrons  by  J.  Jackson  (  op  cit  pp. 
619-637  ): 


dE  _  4rce4 
dx 


me 


NZbC' 


lund 


V 


3 


^  r  , — 

inlfy-il  fiEL 

)  \J  2 


me 

h<co> 


1 

.1 

J 


(in-5) 


where  I  *  h<cu>. 


These  formulae  differ  only  in  the  arguments  of  the  logarithm  and  the  correction  terms 
involving  |3.  A  comparison  of  this  correction  term  follows. 


Table  III-l:  Comparison  of  the  correction  terms  of  equations  (ID-3)  and  (HI-4) 


3 

p2 

-±S 

2  w 

7* 

-iln2+iA 

2  2 

■ 

-  0.005 

mm 

. 

-  0.020 

-H9 

-0.21 

I'Jv 

-  0.045 

0.11 

-0.24 

0.4 

. 

-  0.080 

0.08 

-0.27 

0.5 

0.25 

-  0.125 

0.04 

-0.31 

0.6 

0.36 

-  0.180 

-0.01 

-0.36 

0.7 

0.49 

-  0.245 

-0.06 

-0.41 

O.S 

0.64 

-  0.320 

-0.10 

-0.45 

0.9 

0.81 

-  0.405 

-0.12 

-0.47 

Since  the  value  of  the  logarithm  term  was  of  order  10  for  the  initial  energies  of  the 
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positrons,  I  considered  the  differences  in  these  small  terms  unimportant 


nib)  Plasma  excitation 

The  positrons  also  interact  with  the  free  electrons  in  a  medium.  Bussard  claims  that 
positrons  lose  energy  to  the  ionized  component  of  a  medium  by  exciting  plasma  waves 
(Bussard,  Ramaty,  and  Drachamn  1979).  The  formula  that  he  quotes  from  Book  and  Ali 
(1975)  is : 


dE 

dx 


1.3x10  13 

E 


f  E  ) 

r  \  i 

E 

1  7=7  1  • 

lkT  J 

M' 

kT 

InA 


eV/cm 


(m-6) 


E 

kT 

where  M  =  _2_  f  dx/x  e'x  , 

■IH  J 
o 

ne=  the  free  electron  density  =  NZ^, 
E  =  the  positron  energy  in  eV, 


f  kT  1 

2 

m 

V  k  1 

' 

A  = 

,4™ee2. 

max 

^.u2  •  J 
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Bussard  et  al.  went  into  no  details  about  the  origin  of  this  equation.  Jackson  deals  with 
the  problem  of  an  ion  of  charge  ze  passing  through  an  electronic  plasma  (op  cit  pp. 
641-643).  The  interactions  are  broken  into  two  types  depending  on  whether  the  impact 
parameter  is  larger  than  or  smaller  than  a  Debye  length.  When  the  impact  parameter  is 
smaller  than  a  Debye  length  there  is  a  two-body  screened  potential  interaction.  When  the 
impact  parameter  is  larger  than  a  Debye  length  there  is  a  collective  response  of  the 
medium  and  plasma  oscillations  are  excited  where  the  energy  of  these  oscillations  is 
extracted  from  the  particle.  Thus  Jackson's  formula  is  broken  into  two  parts: 
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After  making  the  appropriate  substitutions: 
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where^^-=  5.1X10"25  MeV-cm2, 
me2 

and  as  before  ne=  the  free  electron  density  =  NZ^. 


In  order  to  make  these  formulae  appropriate  for  positrons  I  believe  that  the  form  of  b^ 
should  be  altered  to  that  of  an  electron: 


V  J 


Jackson  states  that  his  calculations  are  for  non-relativistic  particles.  Bussard  makes  no 
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statement  about  the  energies  where  his  formula  is  valid.  Jackson's  small  impact 
parameter  formula  and  Bussard's  formula  are  nearly  equal  if  one  assumes  that  the  latter  is 
also  non-relativistic.  But  the  problem  is  that  they  seem  to  be  explaining  two  different 
types  of  interactions,  plasma  oscillations  and  two-particle  interactions.  This  is  one  of  my 
major  dilemmas,  considering  that  plasma  energy  loss  is  not  the  objective  of  this  thesis. 
What  follows  is  a  comparison  of  the  coefficients  of  Bussard's  and  Jackson’s  formulae 
(assuming  both  formulae  are  valid  at  the  listed  energies ). 


Table  HI-2: 

Comparison  of  the  coefficients  in  ec 

uations  (HI-6)  and  (UI-7) 

E  (MeV) 

1.3xH) —  Me  V -cn? 

1.50 

ESI 

0.87 

5.5 

1.00 

1.3 

5.8 

0.80 

m 

1.6 

6.0 

0.60 

0.79 

2.2 

6.5 

0.40 

0.69 

3.3 

7.4 

0.20 

0.48 

6.5 

11 

0.10 

0.30 

13 

17 

0.17 

26 

30 

0.04 

130 

130 

25 

(all  values  are  xltf  ) 


The  total  energy  loss  due  to  the  ionized  component  is  either  equation  (HI-6): 
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or  by  adding  equations  (HI-7)  and  (HI- 8): 
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Hie)  Radiative  losses  due  to  Bremsstrahlung 

The  radiative  loss  rate  due  to  Bremsstrahlung  has  been  calculated  by  several  people 
also.  Bethe  and  Heitler  calculated  the  loss  rate  to  be  (Wu  1960): 
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where  me2  «  EQ  «  -  ^  and  screening  of  the  nucleus  is  neglected. 


The  ratio  of  radiative  energy  loss  to  ionization  loss  is  given  approximately  by: 
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=  0.0685xE  (MeV). 
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Thus  this  formula  predicts  radiative  losses  to  be  small  at  our  energies. 

Jackson  also  has  words  of  wisdom  on  this  subject  (op  cit  pp.  708-719).  He  gives 
energy  loss  equations  for  the  nonrelativistic  limit  and  the  ultrarelativistic  limit .  The  ratios 
of  radiative  losses  to  ionization  losses  are: 
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where  B  is  the  argument  of  the  logarithm  in  the  ionization  energy  loss  equation. 


These  formulae  yield  ratios  of  a  few  ten  thousandth's  of  a  percent  and  a  few  percent 
respectively.  Zombeck  also  has  a  formula  for  radiative  losses  (op  cit  Sect.  13  ,  p.  23): 

=  3.09xl0‘27  NZ2  ( E  +  me2 )  Mev/cm  .  (IH-14) 

dx 

This  is  also  negligibly  small  compared  to  the  ionization  loss  rate  for  our  case.  Thus  I  feel 
that  Bremsstrahlung  can  be  safely  ignored  in  our  calculations.  I  believe  that  Cherenkov 
can  also  be  ignored  because  of  the  small  rate  of  loss  (Wu  1960),  if  it  exists  at  all. 

nid)  Conclusion 

From  these  results  I  must  choose  expressions  for  the  energy  loss  of  56Co  positrons 
in  the  supernova  interior.  Because  the  fraction  of  positions  that  survive  the  expansion  is 
only  weakly  sensitive  to  the  differences  in  thermalization  times  the  choice  of  which  to  use 
is  not  crucial.  I  chose  to  use  equation  (ID-4)  and  equation  (IH-9)  (see  figures  4  and  5). 


Ionization 


(uio/asW)  2 1  ©I  x  xp/ap 


A  plot  of  the  ionization  loss  (equation  (III-4))  wuh  the  number  of  bound  electrons  per  atom 
set  at  24.5  and  the  density  of  atoms  set  at  1 .08x  1 0  atoms/cm  .  Note  the  large  increase  in  dE/dx 

as  the  energy  falls  below  0.2  MeV. 


IV.  Hubble  Flow  and  Its  Impact  on  the  Density  Evolution  of  Supernovae 


As  I  have  shown,  the  energy  loss  rate  for  a  positron  passing  through  matter  depends 
on  the  density  of  that  matter.  The  rate  of  positron  annihilation  also  depends  on  the 
density  of  the  matter,  specifically  the  electron  density.  Thus  I  had  to  find  an  expression 
for  the  density  evolution  of  the  supernova  ejecta.  An  important  step  in  simplifying  the 
density  evolution  of  a  spherical  shell  in  a  supernova  was  the  assumption  that  the 
expanding  gas  is  in  a  Hubble-type  flow.  Hubble  flow  simply  means  that  the  velocity  at  a 
given  point  in  the  supernova  is  proportional  to  its  distance  from  the  center  of  the 
explosion.  It  is  easy  to  understand  why  this  type  of  outward  flow  should  occur.  As  the 
shock  propagates  outward  each  shell  is  pushing  on  the  layers  on  top  of  it.  The  outer 
shells  are  pushing  on  less  matter  than  inner  ones  and  therefore  should  accelerate  to  a 
larger  velocity.  This  becomes  true  in  supemovae  after  acceleration  has  ended  and  after 
the  size  is  much  greater  than  the  initial  size.  From  this  time  on  the  ratio  of  velocity,  v(r), 
to  distance  from  the  center,  r,  is  a  constant  for  all  r.  This  also  means  that  there  is  a 
velocity  gradient  across  a  spherical  shell.  If  one  considers  such  a  shell,  it  can  be  shown 
as  follows  that  the  density  within  that  shell  decreases  as  r3. 
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As  t  =>  t  +  dt:  r  =>  r  +  vdt 

r  +  dr  =>  r  +  dr  +  (v  +  dv)dt. 

From  this,  one  can  show  that  the  thickness  of  the  shell  has  increased  to  dr  +  dvdt  =  dr(l  + 
(VR/R)dt)  where  VR  is  the  velocity  of  the  outermost  region  of  the  supernova  and  R  is  the 
radius  of  that  region.  It  is  easy  to  show  dv  =  dr  (Vr/R)  because  in  a  Hubble  flow  v(r)/r  = 
Vr/R  =  constant.  To  show  that  the  density  in  this  region  is  proportional  to  r3  one  needs 
to  use  conservation  of  mass  within  the  shell: 

M(t)  =  M(t  +  c!t), 

where  M(t)  =  4rcr(t)2  dr(t)  p(t),  and 

M(t  +  dt)  =  4ttr(t  +  dt)2  dr(t  +  dt)  p(t  +  dt). 

Setting  M(t)  =  M(t  +  dt): 

==>  4rc  r(t)2  dr(t)  p(t)  =  4tc  r(t  +  dt)2  dr(t  +  dt)  p(t  +  dt). 

Keeping  terms  only  to  first  order  in  dt: 

47tr(t)2dr(t)  p(t)  =  47t  [r(t)2  +  2r(t)v  dt]  dr(t)  [1  +  (Vr/R)  dt]  p(t  +  dt), 

=>  r(t)2  p(t)  =  [r(t)2  +  2r(t)v  dt]  [1  +  (Vr/R)  dt]  p(t  +  dt), 

==>  0  =  r2  dp  +  [2rv  +  r2(VR/R)]  p  dt , 

=  r2  [dp  +  3(Vr/R)  p  dt], 

=>  dp  =  -  3(Vr/R)  p  dt 

The  solution  to  this  equation  is: 

t30 

P(t)  =  p(to)  -f ,  where  t  >  ^ , 
t 
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or 

4 

n(t)  =  n(tQ)  ,  where  t>tQ.  (IV- 1) 

Thus  the  velocity  gradient  across  a  spherical  shell  due  to  a  Hubble-type  flow  causes  the 
particle  density  within  that  flow  to  decrease  as  r3  with  the  restriction  that  the  value  taken 
for  to  must  be  after  a  Hubble-type  flow  has  been  established.  For  my  study  the  values 
n(to)  and  time  to  are  obtained  from  published  supernova  models  (see  section  VI). 
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1L _ The  Radiative  Capture  and  Annihilation  of  a  Positron 

So  far  I  have  covered  the  production  and  thermalization  of  positrons.  Positrons  are 
destroyed  by  annihilating  with  electrons.  Positrons  can  annihilate  directly  or  by  forming 
positronium  first.  Positronium  in  the  singlet  state  has  a  lifetime  of  10*10  s  and  emits  two 
511  keV  photons  during  annihilation.  Triplet  positronium  has  a  lifetime  of  10-7s  and 
undergoes  three-photon  annihilation  (Bussard,  Ramaty,  and  Drachman  1979;  Ramaty  and 
Lingenfelter  1987).  The  probability  of  annihilation  in  flight  is  negligible  until  the 
positrons  slow  to  energies  of  several  hundred  eV  (Bussard,  Ramaty,  and  Drachman 
1979).  The  values  of  equations  (HI-4)  and  (HI-9)  rise  dramatically  at  low  energies 
(below  0.2  MeV).  The  positrons  will  therefore  spend  only  a  small  amount  of  time  at 
these  energies  before  thermalizing.  I  have  therefore  required  all  positrons  to  theimalize 
before  annihilating.  At  temperatures  associated  with  the  interior  of  a  supernova  (a  few 
thousand  degrees)  the  dominant  mode  of  annihilation  is  positronium  formation  by 
radiatively  recombining  with  free  electrons  (Bussard,  Ramaty,  and  Drachman  1979)  (see 
figure  6).  The  cross  section  for  such  a  reaction  changes  with  velocity  (and  therefore 
temperature).  Thankfully  the  cross  section  is  only  proportional  to  T_1/2  (Ostcrbrock 
1974)  and  is  therefore  relatively  insensitive  to  these  changes. 

The  theory  of  radioactive  decay  is  well  known.  The  radiative  recombination  and 
subsequent  annihilation  of  a  positron  can  be  treated  in  an  analgous  manner.  The 
probability  that  a  positron  will  decay  (annihilate)  in  the  next  time  interval  dt  is 
proportional  to  dt  and  inversely  proportional  to  the  mean  lifetime  t : 

Probability  of  Decay  *  dt/x  .  (V-l) 


If  the  mean  lifetime  is  taken  to  be  a  constant,  i.e.  constant  density  and  constant 
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recombination  rate,  and  if  no  new  positrons  are  created  then  the  number  of  positrons 
decreases  exponentially  with  time: 


N+  =  N+0e'x.  (V-2) 

In  an  expanding  supernova  the  density  is  proportional  to  t*3,  as  I  have  shown.  Since  the 
mean  lifetime  against  positronium  formation  is  inversely  proportional  to  the  density,  the 
mean  lifetime  goes  as  t3  : 


T 


=  l 20  ^free 


>1 
X  I 
J 


(V-3) 


where  to  =  a  time  when  Hubble  flow  has  been  established, 
n(to)  =  the  density  of  atoms  at  to, 

=  the  number  of  free  elections  per  atom  at  to, 

\  =  the  radiative  recombination  rate, 

jvr 

therefore  -rrt-=--^T  •  (V-4) 

N+  at3 


There  are  then  two  ways  of  calculating  the  probability  of  a  positron  surviving  to  a 
predetermined  time: 

1)  A  Monte  Carlo  simulation  using  the  probability  of  capture  and  annihilation, 

2)  An  analytic  solution  to  equation  (V-4)  which  requires  taking  a  as  a  constant 


26 


The  Monte  Carlo  routine  that  I  used  stepped  through  time  with  a  constant  time  step. 
At  the  beginning  of  each  step  the  mean  lifetime  was  evaluated  for  that  time.  This  set  the 
probability  of  capture  for  that  time  step.  A  random  number  was  generated  for  each 
positron  that  had  survived  to  that  time.  If  the  number  (between  0  and  1)  was  less  than  the 
probability,  the  positron  was  removed.  The  time  was  then  incremented  after  all  the 
surviving  positrons  had  been  checked  for  survival.  Thus  the  number  of  positrons  that 
survived  to  the  "end"  time  out  of  the  initial  number  injected  could  be  found.  There  was 
one  major  difficulty  with  this  method.  The  "resolution"  was  inversely  proportional  to  the 
number  of  positrons  initially  injected.  By  "resolution"  I  mean  the  smallest  probability  of 
survival  obtainable.  For  example,  in  following  the  history  of  104  particles  the  probability 
of  survival  cannot  be  calculated  when  it  becomes  comparable  to  1(H  or  less.  The 
survival  fraction  is  quantized  in  units  of  104.  In  Type  II SN  the  probability  of  survival  is 
a  very  small  number.  To  get  probabilities  down  to  the  size  needed  would  use  an 
enormous  amount  of  computer  time.  This  problem  led  to  the  use  of  the  following 
analytic  solution. 

The  analytic  solution  I  derived  from  equation  (V-4)  proved  to  be  more  useful  for  my 
calculations : 


(V- 5) 
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With  this  formula  I  could  calculate  the  probability  of  a  positron  surviving  without  using 
arbitrarily  large  numbers  of  initial  positrons  to  get  small  probabilities. 

A  comparison  of  the  two  methods  shows  that  they  produce  comparable  results  when 
the  number  of  positrons  that  survive  to  the  end  of  the  monte  carlo  program  is  of  the  order 
100  (see  Appendix  A).  Equation  (V-5)  can  be  used  to  produce  an  analytic  model  for  the 
survival  of  positrons  but  it  is  difficult  to  incorporate  the  thermalization  process  (see 
Analytic  Models  of  Positron  Survival,  section  VII). 

The  advantage  of  the  numerical  approach,  on  the  other  hand,  is  that  it  does  not 
require  a  constant  annihilation  rate  per  free  electron.  The  temperature,  ionization,  direct 
in-flight  annihilation,  etc.  can  all  be  represented  by  time-dependent  quantities  in  a  Monte 
Carlo  history.  That  approach  is  therefore  more  general.  The  temperature  of  the 
supernova  ejecta  should  decrease  with  time  due  to  the  fewer  number  of  decays  depositing 
energy  within  the  supernova.  The  rate  of  this  cooling  is  unknown  because  it  is  uncertain 
whether  a  pulsar  was  created  in  the  SN  1987A.  This  and  other  unknowns  (degree  of 
radial  mixing,  expansion  velocity,  magnetic  field  configuration,  clumping,  etc.)  coupled 
with  the  weak  temperature  dependence  of  the  recombination  rate,  X,  led  me  to  believe  that 
my  analytic  solution  was  sufficient  for  this  rough  survey  of  positron  survival.  I  have 
therefore  used  the  approximation  that  the  mean  lifetime  of  positrons  against  radiative 
recombination  is  exactly  proportional  to  t3  and  have  esimated  the  proportionality 
constants  for  typical  supernova  conditions  in  my  models. 
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YL _ Models  of  SN  1987 A 

I  chose  to  use  Type  II  supernova  models  by  Woosley  et  al.  and  Nomoto  et  al.  as  the 
sources  of  parameters  needed  as  input  in  my  models.  The  model  of  Woosley’s  that 
produced  the  best  fit  to  the  light  curve  of  SN  1987A  was  model  10HM.  This  model 
assumed  that  the  progenitor  star  was  Sk  -69  202.  This  star  was  assumed  to  have  been  20 
M0  when  on  the  main  sequence.  Such  a  star  would  have  a  6  M0  helium  core.  The  star 
was  then  allowed  to  lose  4  M0  from  its  hydrogen  shell  prior  to  the  explosion.  A  1.4  M0 
neutron  star  was  formed  from  this  helium  core.  The  explosion  was  simulated  by 
replacing  the  neutron  star  with  a  piston  that  had  an  energy  of  1.45x10s1  ergs.  During  the 
explosion,  approximately  0.07  M0  of  56Ni  was  formed  in  the  silicon  shell.  This  model 
mixed  the  56Ni  radially  outward  to  a  mass  coordinate  of  9  M0.  This  was  done  to 
improve  the  light  curve  fit  in  the  20  -  40  day  region.  It  also  accounted  for  the  early 
detection  of  x-  and  y-rays  at  175  days.  The  regions  where  the  56Ni  was  mixed  was  not 
homogenized  and  the  model  lOH(non-mixed)  density  profile  was  left  unaltered  (Woosley 
1988;  Woosley  (in  press);  Pinto  and  Woosley  1988). 

Nomoto  et  al.  have  also  produced  several  hydrodynamic  models  of  SN  1987 A.  The 
model  that  produced  the  best  fit  to  the  observed  light  curve  was  model  1 1E1.  This  model 
also  assumed  that  the  progenitor  star  was  20  M0  while  on  the  main  sequence.  Similarly, 
a  1.4  M0  neutron  star  was  formed  from  the  6  M0  helium  core.  In  this  model  the  star  was 
assumed  to  have  undergone  mass  loss  leaving  it  with  12.7  M0  of  material  at  the  time  of 
collapse.  The  explosion  was  simulated  by  replacing  the  iron  core  with  a  point  mass  and 
introducing  the  explosion  energy  interior  to  the  1 1.3  M0  of  ejecta.  Model  1 1E1  used  an 
explosion  energy  of  l.OOxlO51  ergs.  This  model  produced  appproximately  0.07  M0  of 
S6Ni.  Radial  mixing  was  also  required  in  this  model  to  fit  observations  of  the  light  curve 
between  25  -  40  days  and  the  early  detection  of  x-  and  y-rays.  Thus  Nomoto  mixes  the 


30 

nickel  out  to  a  radial  mass  coordinate  of  10  M0  (Note:  Nomoto  and  Woosley  use  different 
conventions.  Nomoto 's  mass  coordinate  system  has  its  zero  point  at  the  inner  edge  of  the 
ejecta  while  Woosley's  is  at  the  center  of  the  remnant.  Hence  forth,  I  will  always  add  1.6 
M0  onto  Nomoto's  coordinate  for  ease  in  comparing  results.)  (Nomoto  et  al.  1988; 
Shigeyama,  Nomoto,  and  Hashimoto  1988). 


Table  VI-1:  Summary  of  supemovae  model  parameters 


Nomoto  et  al. 

Woosley  et  al. 

Main  Sequence  Mass 

20  M0 

20  M0 

Neutron  Star  Mass 

1.4  MQ 

1.4  Mg 

Mass  of  Ejecta 

11.3  M0 

14.6  Mq 

Energy  of  Explosion 

l.OOxlO51  ergs 

1.45xl051  ergs 

Maximum  Radius 

for  Mixing  Ni 

11.6  M0 

(Woosley  Conv.) 

9Mg 

The  quantities  I  required  from  these  models  were  a  density  profile  at  an  early  time 
after  the  explosion  (but  after  Hubble  flow  had  been  established),  a  mixing  profile,  and  an 
approximate  temperature  of  the  region.  Fortunately  all  the  needed  information  was 
available  (see  figures  7-10).  For  each  model  I  divided  the  region  of  the  ejecta  which 
contained  nickel  into  shells.  I  approximated  the  density  profiles  for  each  model  by  a 
series  of  step  functions  to  obtain  a  typical  density  for  each  shell.  I  determined  an  analytic 
function  that  fit  Woosley's  mixing  profile  and  integrated  it  over  the  region  of  each  shell  to 
determine  the  amount  of  nickel  in  that  shell.  Nomoto's  mixing  was  a  step  function  so  the 
amount  of  nickel  in  each  of  those  shells  was  just  the  area  under  each  step  function.  A 
summary  of  the  results  follows. 


Table  VI-2:  Summary  of  model  10HM  (Woosley  et  al.) 


Shell  (M0) 

56  Ni  (M0 ) 

n(tg)  (cm)3 

p(t0)  (g-cm)3 

1.6  -  2.0 

0.0105 

1.08xl010 

l.OxlO*12 

2.0  -  2.5 

0.0109 

5.38xl010 

5.0x10  42 

2.5  -  3.0 

0.0090 

3.76xl0n 

3.5xl041 

3.0  -  3.5 

0.0073 

5.38xl010 

5.0x10  42 

3.5  -  4.0 

0.0061 

3.76x10 10 

3.5x10  42 

4.0  -  4.5 

0.0050 

2.37xl010 

2.2xl042 

4.5  -  5.0 

0.0041 

2.37xl010 

2.2xl042 

5.0  -  5.5 

0.0034 

2.16xl010 

2.0xl042 

5.5  -  6.0 

0.0027 

1.61xl010 

1.5xl042 

6.0  -  6.5 

0.0022 

3.23xl09 

3.0x10  43 

6.5  -  7.0 

0.0019 

2.16xl09 

2.0x10  43 

7.0  -  7.5 

0.0015 

2.69x1 09 

2.5xl043 

7.5  -  8.0 

0.0012 

3.23xl09 

3.0x10  43 

8.0  -  8.5 

0.0009 

3.23xl09 

3.0x10  43 

8.5  -  9.0 

0.0008 

3.23xl09 

3.0x10  43 

(tjj  =  86  days) 


Table  VI-3:  Summary  of  model  1 1E1  (Nomoto  et  al.) 


Shell  (Mgj ) 

56  Ni  (Ms  ) 

n(tn)  (cn?J3 

p(tn)  (g-cmJ3 

1.6  -  2.1 

0.0100 

1.08xl010 

l.OxlO42 

2.1  -  2.6 

0.0079 

8.60xl010 

8.0x10  42 

2.6  -  4.0 

0.0222 

1.08xl010 

l.OxlO42 

4.0  -  5.1 

0.0155 

9.68xl09 

9.0x10  43 

5.1  -  6.2 

0.0078 

5.38xl09 

5.0x10  43 

6.2  -  7.6 

0.0056 

9.68x10* 

9.0x10  44 

7.6  -  9.6 

0.0008 

4.30x10* 

4.0x10  44 

9.6-11.6 

0.0007 

2.15x10* 

2.0x10  44 

(where  tg  =  116  days) 
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Figure  7:  The  density  profile  of  Woosley’s  model  10H  at  various  times  after  the  explosion.  I 
used  die  profile  at  86  days  for  my  initial  conditions  (Woosley  1988). 


The  mixing  profile  for  Woosley's  model  10HM.  This  profile 
amount  of  cobalt  in  each  shell  (Woosley  (in  press)). 


34 


(£**W3/9) OHy 


M  (R) /MSUN 

Figure  9:  The  density  profiles  for  Nomoto's  models  1 1E1  (solid)  and  1 1E2  (dashed)  at 
t  *  1 16  days  (Shigeyama,  Nomoto,  Hashimoto  1988). 


The  mixing  profile  for  Nomoto’s  model  1 1E1.  This  profile  was  used  to  calculate  the 
amount  of  cobalt  in  each  shell  (Nomoto  et  al.  1988). 
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VTI.  Analytic  Models  of  Positron  Survival 

I  have  worked  on  an  analytic  model  that  began  with  a  differential  equation  involving 
the  time  rate  of  change  of  the  number  of  positrons  per  gram.  It  contained  a  source  term 
due  to  the  production  of  positrons  by  the  decay  of  56Co  and  a  loss  term  due  to  the 
radiative  recombination  of  positrons  with  free  electrons.  The  main  problem  with  this 
solution  is  that  it  is  hard  to  incorporate  the  process  of  thermalization.  But  I  have  found 
that  the  time  it  takes  for  a  positron  to  become  thermal  is  small  cc  npared  to  time  of 
injection  into  the  supernova  (at  least  for  injection  times  of  interest).  Therefore  the  error 
from  this  shortcoming  should  be  small. 

The  analytic  solution  begins  with  the  expression: 

^iUnen+=ScV.  (VII-l) 
where  n+  =  the  number  of  positrons  per  gram, 

S  =  initial  56Co  decay  rate  per  gram  =  (NAvog  /56)  x  X^, 

Xq  =  56Co  decay  constant  =  1/111.52  days, 

ne  =  the  number  of  free  electrons  per  cm3  =  the  number  of  atoms  per  cm3  x  Z^, 
X  =  the  radiative  recombination  rate. 

But  from  equation  (IV- 1)  the  density  of  atoms  is: 

lo 

n(t)  =  n(t0)  -j . 

This  changes  equation  (VII-l)  to: 


+  (X  nCtjj)  )  t 3  n+  *  S  e  ^ 


(vn-2) 
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An  equation  of  the  form  dy/dx+a(x)y  =  h(x)  is  equivalent  to  the  equation  d(py)/dx=ph  if : 


p  =  exp  J  a(x)  d*  . 


The  solution  of  this  equation  is  (Hildebrand  1976): 


y  =  j  Jphdx  + 


const. 

P 


In  particular,  if  y  =  yi  at  x  =  xx: 


yW  =  jhW£|ldx'  +  y,^.  (vn-3) 


By  making  the  substitution  in  equation  (V-3): 


aC)=Xn(t0)^Zfreet-3=-iJ, 

at 

r 

■>  p(0 *cxP^*2a(t2  Jj • 


At  t  =  tlf  n+(t)  =  n+(tt)  for  tL  ^  (the  initial  time  of  the  expansion).  From  equation 


(VH-3): 


"  i  fx-xT 

n+(t)  =  n+(tl)exp|^2^^t2  t\  + 


r  rnl  t  r  f  . 

<pl  t2  ]  I  [  expl  -  r°  1  +  f)2  I  dt' . 

KL2aUJjJ  \  l  2at  JJ 


(vn-4) 
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The  first  term  corresponds  to  the  number  of  positrons  existing  at  tt  that  survive  until  time 
t.  The  second  term  corresponds  to  the  number  of  positrons  emitted  between  tt  and  t  that 
survive  until  time  t.  I  then  let  tj  =>  to  and  let  t  =>  °o  ; 


r  i  i 


n+(°°)  =  n+(t0)  exp  I  "  2at2Q  j  +  S  J  exp  I  - 


L  l 


V  + 


i 


'll 


2at’2JJdt' 


(vn-5) 


The  integral  needed  to  be  done  numerically.  To  do  this  I  used  a  routine  on  the  microVax, 
DQDAGI,  from  the  IMSL  Math  Library.  When  time  was  kept  in  seconds  it  was  difficult 
to  get  the  integrating  routine  to  produce  a  non-zero  solution  to  the  integral  in  equation 
VH-5  because  of  the  extremely  small  size  of  the  integrand  over  the  entire  range  of  the 
integration  except  for  a  very  sharp  peak  at  5xl08  seconds  far  from  the  lower  limit  of 
integration.  I  therefore  converted  my  units  of  time  into  years.  This  moved  the  peak 
closer  to  the  lower  limit  of  integration  (in  the  absolute  sense)  and  thus  made  it  harder  for 
the  integration  routine  to  ignore.  I  used  this  model  to  calculate  the  number  of  positrons 
that  would  survive  until  t  =  «o  from  the  innermost  shell  of  Woosley's  model  10HM.  I 
used  the  following  values  for  the  input  parameters: 


Xo  =  1.038xl0-7  s-1  =  1/111.52  days  =  3.27  yr-i, 

X  =  2.00x10-12  cm3  s-i  =  6.31xl0’5  (T-3500  K)  (Bussaid  et  al  1979), 
Zfree=  1.5  electrons/atom, 
to  =  7.43xl06  s  =  86  days  =  0.236  yrs, 
n(to)  =  1.08xl010  atoms  cm-3, 

a  =  7.52x10-2°  s-2  -  7.45xl0*5  yrs-2. 


Using  these  parameters,  the  value  for  the  integral  in  equation  VH-5  was  8.52xl0*27  s. 
The  first  term  in  equation  (VH-5)  because  it  was  very  small.  By  taking  a  value  for  S  of 
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1.12xl015  g-1  s'1, 1  calculated  that  9.54xl042  positrons  per  gram  should  survive.  This 
translates  to  3.99xl019  positrons  surviving  from  the  inner-most  shell  of  Woosley's 
model  10HM.  This  result  is  compared  to  that  of  the  computer  simulation  in  section  IX. 
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Vni.  Computer  Model  of  Positron  Survival 

My  computer  model  for  positron  escape  from  a  Type  II  supernova  drew  from  the 
derivations  and  research  presented  in  sections  II  -  VI  of  this  thesis.  Two  of  the  major 
improvements  of  this  type  of  model  over  the  previously  discussed  analytic  model  were: 

1)  Requiring  positrons  to  thermalize  before  annihilating, 

2)  Allowing  the  radial  mixing  of  cobalt  into  regions  of  lower  density. 

To  accomodate  the  radial  mixing,  the  supernova  ejecta  was  broken  up  into  a  number  of 
shells  (as  discussed  in  section  VI).  There  were  two  sets  of  initial  conditions 
corresponding  to  the  two  supernova  models  discussed  in  section  VI.  Each  shell  was 
assigned  an  initial  density  and  amount  of  56Ni  (56Co)  in  accordance  with  the  model  being 
used.  A  summary  of  this  data  was  shown  in  tables  VI-2  and  VI-3.  Other  parameters 
needed  in  this  model  were  the  56Co  decay  rate,  the  radiative  recombination  rate  for 
positrons  and  electrons,  and  the  degree  of  ionization  of  the  ejecta.  The  values  I  took  for 
these  parameters  were  the  same  as  I  used  in  the  analytic  model: 


Xq=  1.038x10-7  s*1, 

2.00xl0'12  cm3  s-i, 

Zfree  =1.5  electrons/atom. 

This  program  calculated  the  number  of  positrons  that  survived  from  each  shell.  To 
do  this,  the  continuous  emission  of  positrons  was  approximated  by  a  discrete  set  of 
positron  injection  times  varying  from  5-30  years  with  the  interval  between  injection 
times  being  1  year.  The  number  of  cobalt  decays  that  occurred  within  each  time  interval 
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dt  centered  on  the  injection  time  can  be  described  by  the  expression: 

dNCo  =  NQ  Co  J^exp  ('  V  1 '  )  j .  exp  ('  ( 1  +  ^t/2  ) )  j  . 

The  number  of  positrons  injected  was  1/5  this  number. 

For  each  injection  time  the  program  calculated  a  mean  time  of  thermalization  for  a 
positron  (the  time  when  a  typical  positron  was  slowed  to  thermal  energies).  This  was 
done  by  using  equations  (HI-4)  and  (ID-9)  to  calculate  a  time  of  thermalization  for  a 
positron  in  each  of  the  energy  intervals  listed  in  table  II- 1.  These  times  were  then 
weighted  by  the  energy  distribution  function  to  get  the  mean  time  of  thermalization. 
During  the  thermalization  process  the  density  was  allowed  to  evolve  according  to 
equation  (IV- 1).  The  program  then  used  equation  (V-5)  to  calculate  the  fraction  of 
positrons  surviving  from  the  mean  time  of  thermalization  until  t  =  °°  and  that  time  when 
the  mean  lifetime  of  positrons  to  radiative  recombination,  t,  equalled  100  years  (as 
calculated  for  the  central  shell  of  model  10HM).  This  second  calculation  was  made  to 
estimate  the  error  in  terminating  the  recombination  process  at  approximately  110  years 
(when  t  =  100  years).  This  could  be  of  use  when  performing  a  Monte  Carlo  history. 
Thus  the  number  of  positrons  surviving  from  a  given  shell  with  a  given  injection  time 
could  be  calculated  by  taking  the  product  of  the  survival  fraction  and  the  number  of 
positrons  injected  for  that  injection  time  in  that  shell.  The  program  stepped  through 
injection  times  of  5  to  30  years  and  through  all  of  the  shells.  The  program  added  the 
contributions  from  all  times  and  all  shells  to  get  a  total  number  of  positrons  surviving  for 
each  model.  A  flow  diagram  for  the  program  along  with  a  printout  are  included  in 
Appendix  B.  A  quantity  of  potentially  more  interest  was  the  probability  of  a  positron 
surviving  from  a  particular  shell.  This  gives  insight  into  how  the  survival  of  positrons 
depends  on  the  density  of  a  region. 

A  summary  of  the  results  for  each  shell  when  the  positrons  were  required  to  survive 
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until  t  =  °°  can  be  seen  in  tables  VIII- 1  and  VICI-2.  Note  the  large  difference  from  shell  to 
shell  in  the  probability  of  surviving.  This  is  due  to  the  fact  that  the  probability  of 
surviving  is  an  exponential  involving  the  density  and  one  can  see  that  the  density  varies 
significantly  from  shell  to  shell.  Thus  the  radial  mixing  of  nickel  (cobalt)  outward 
significantly  improves  the  production  of  positrons  able  to  contribute  to  a  diffuse 
annihilation  line  (see  figures  11  and  12). 

As  I  have  already  stated  I  have  also  calculated  the  probability  that  a  positron  will 
survive  until  the  mean  lifetime,  x,  equalled  100  years  (approximately  110  years  after  the 
explosion)  and  the  number  that  will  do  so  from  each  shell  to  compare  with  the  previous 
results  (see  tables  VTII-3  and  Vm-4).  The  differences  between  the  results  for  survival 
until  t  =  oo  and  until  x  =  100  years  is  small  for  the  shells  producing  the  most  positrons. 
The  difference  between  the  overall  survival  probability  and  the  overall  number  of 
"survivors"  is  correspondingly  small,  1 1.4%  for  model  10HM  and  2.6%  for  model  1 1E1 
(see  table  VIII-5). 


Table  VIII-1:  Summary  of  Results  for  model  10HM  Initial  Conditions 


Shell  (Mo ) 

Initial  Density 
(at  86  days) 

Probability  of  Surviving 
until  t  =  °o 

Number  That  Survive 

until  t  =  oo 

1.6  -  2.0 

1.08xl010 

1.98xl0'33 

8.89xl019 

2.0  -  2.5 

5.38xl010 

1.87xl0'57 

0 

2.5  -  3.0 

3.76xlOn 

9.52xl0'155 

0 

3.0  -  3.5 

5.38xl010 

1.87xl0-57 

0 

3.5  -  4.0 

3.76xl010 

6.53xl0'51 

1.71xl02 

4.0  -  4.5 

2.37xl010 

1.40x10- 43 

3.00x1 09 

4.5  -  5.0 

2.37xl010 

1.40xl043 

2.46xl09 

5.0  -5.5 

2.16xl010 

3. 08x10  42 

4.48xl010 

5.5  -  6.0 

1.61xl010 

2.98xl0‘38 

3.44xl014 

6.0  -  6.5 

3.23xl09 

3.16xl0'22 

2.97xl030 

6.5  -  7.0 

2.16xl09 

2.14xl0'19 

1.74xl033 

7.0  -  7.5 

2.69xl09 

6.83xl0'21 

4.38xl031 

7.5  -  8.0 

3.23xl09 

3.16xl0*22 

1.62xl030 

8.0  -  8.5 

3.23xl09 

3.16xl0'22 

1.22xl030 

8.5  -  9.0 

3.23xl09 

3.16xl0'22 

1.08xl030 

Table  VIII-2:  Summary  of  Results  for  model  1 1E1  Initial  Conditions 


Shell  <M  ) 

Initial  Density 
(at  116  days) 

Probability  of  Surviving 
until  t  =  <» 

Number  that  Survive 

until  t  =  oo 

1.6  -  2.1 

1.08xl010 

3.04xl045 

1.30x10s 

2.1  -  2.6 

8.60x10 10 

1.80xl0‘105 

0 

2.6  -  4.0 

1.08xl010 

3. 04x10 45 

2.88x10s 

4.0  -  5.1 

9.68xl09 

1.32xl043 

8.76xl09 

5.1  -  6.2 

5.38xl09 

9.08xl0'36 

3.03xl017 

6.2  -  7.6 

9.68x10s 

5.01X10'20 

1.20xl033 

7.6  -  9.6 

4.30x10s 

3.27xl0'15 

1.12xl037 

9.6-  11.6 

2.15x10s 

5.06xl0'12 

1.51xl040 
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Table  VIII-3:  Summary  of  Results  for  model  10HM  Initial  Conditions 


Shell  (Mo ) 

Initial  Density 
(at  86  days) 

Probability  of  Surviving 
until  x  =  100  years 

Number  That  Survive 

until  x  =  100  years 

1.6  -  2.0 

1.08xl010 

3.43xl0'33 

1.54xl020 

2.0  -  2.5 

5.38xl010 

2.90xl0'56 

0 

2.5  -  3.0 

3.76xlOn 

2.02xl0‘146 

0 

3.0  -  3.5 

5.38xl010 

2.90xl0’56 

0 

3.5  -  4.0 

3.76xl010 

4.45xlO'50 

1.16xl03 

4.0  -  4.5 

2.37xl010 

4.70x10  43 

1.00x10 10 

4.5  -  5.0 

2.37xl010 

4.70xl043 

8.24xl09 

5.0  -5.5 

2.16xl010 

9.27xl042 

1.35xlOn 

5.5  -  6.0 

1.61xl010 

6.77xl0‘38 

7.82xl014 

6.0  -  6.5 

3.23xl09 

3.72xl0-22 

3.50xl030 

6.5  -  7.0 

2.16xl09 

2.39xl0‘19 

1.94xl033 

7.0  -  7.5 

2.69xl09 

7  83x10 ‘21 

5.03xl031 

7.5  -  8.0 

3.23xl09 

3.72xl0'22 

1.91xl030 

8.0  -  8.5 

3.23xl09 

3.72xl0*22 

1.43xl030 

8.5  -  9.0 

3.23xl09 

3.72xl0*22 

1.27xl030 

Table  VIII-4:  Summary  of  Results  for  model  11E 1  Initial  Conditions 


Shell  (Mq) 

Initial  Density 
(at  116  days) 

Probability  of  Surviving 
until  x  =  100  years 

Number  that  Survive 

until  x  =  100  years 

1.6  -  2.1 

1.08xl010 

l.HxlO44 

5.01x10® 

2.1  -  2.6 

8.60xl0l° 

8.47x10  401 

0 

2.6  -  4.0 

1.08xl010 

l.HxlO44 

l.llxlO9 

4.0  -  5.1 

9.68xl09 

4.44x10 43 

2.94x10 10 

5.1  -  6.2 

5.38xl09 

1.78xl0*35 

5.94xl017 

6.2  -  7.6 

9.68x10s 

5.66X10’20 

1.36xl033 

7.6  -  9.6 

4.30x10® 

3.47xl0’15 

1.19xl037 

9.6-  11.6 

2.15x10® 

5.18xl0'12 

1.56X1040 
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Table  VIII-5:  Summary  of  Overall  Results  for  models  10HM  and  11E1 


10HM 

t  =  oo 

10HM 

x=  100 

11E1 

t  =  oo 

11E1 

t=  100 

Probability  of 
Surviving 

6.20x1  O'21 

6.93xl0'21 

5.03xl0'14 

5.16xl0"14 

Number  of 

Survivors 

1.79xl033 

2.00xl033 

1.52xl040 

1.56xl040 

%  Difference 

in  number  of 

11.7% 

2.6% 

"survivors" 

I  comment  on  these  results  and  compare  them  to  the  analytic  model  in  section  DC. 


This  plot  shows  the  probability  of  a  positron  surviving  until  t  =  °°  with  initial 
conditions  set  from  each  of  the  models.  Note  the  extremely  large  range  of  values 
due  to  variations  in  density. 
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Figure  12:  This  plot  shows  fraction  of  the  total  number  of  surviving  positrons 
contributed  by  each  shell  for  both  sets  of  initial  conditions.  Note  that 
vast  majority  of  the  positrons  are  coming  from  the  shell  with  the  lowest 
density. 
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IX.  Summary  and  Conclusions 

My  original  purpose  for  undertaking  this  project  was  to  investigate  Type  II 
supemovae  as  a  possible  source  for  the  diffuse  component  of  the  511  keV  annihilation 
radiation  from  the  Galactic  center.  At  the  recent  14th  Texas  Symposium  on  Relativistic 
Astrophysics,  Reuven  Ramaty  argued  that  a  source  of  3X1043  positrons/sec  was  needed 
to  account  for  the  diffuse  background  (this  is  easily  verified  to  within  an  order  of 
magnitude  by  a  simple  calculation  if  the  distance  to  the  Galactic  center  is  taken  to  be  8 
kpc).  My  models  show  that  Type  II  supemovae  do  not  produce  nearly  enough  positrons 
to  account  for  this  line.  Using  initial  conditions  set  by  model  11E1,  my  models  show 
only  1.5xl040  positrons  will  survive  (model  10HM  inital  conditions  yielded 
approximately  1033).  If  one  assumes  a  fairly  optimistic  rate  of  one  Type  II  supernova 
every  100  years  in  the  Galactic  center,  a  production  rate  of  only  4.777X1030  positrons/sec 
is  achieved.  This  is  well  short  of  what  is  needed.  This  does  not  mean  that  this  model  was 
a  failure.  There  was  good  agreement  between  the  analytic  model  and  the  computer 
model.  The  computer  model  predicts  8.89xl019  positrons  surviving  until  t  =  °°  from  the 
inner  shell  of  model  10HM  while  the  analytic  model  predicts  3.99xl019.  The  diffence 
can  be  accounted  for  by  the  fact  that  in  the  computer  model  the  positrons  were  required  to 
thermalize  before  annihilating.  During  the  thcrmalization  process  the  density  in  each  shell 
decreased.  This  coupled  with  the  strong  dependence  of  the  probability  of  survival  on  the 
density  increases  the  chance  for  survival  from  a  particular  injection  time.  As  mentioned 
previously,  the  integral  contained  in  the  analytic  solution  was  difficult  to  evaluate.  A 
non-zero  solution  came  only  after  converting  seconds  to  years. 

The  computer  model  produced  interesting  results  as  I  have  stated  in  section  VIII.  The 
most  interesting  is  the  strong  dependence  of  the  probability  of  survival  on  the  density  (see 
figure  13).  The  reason  for  this  is  that  the  probability  of  survival  contains  an  exponential 
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involving  the  density.  Thus  a  variation  in  density  from  shell  to  shell  of  a  factor  of  50, 
produces  a  factor  of  e50  difference  in  the  probability  of  surviving.  This  is  the  reason  that 
mixing  has  such  an  overwhelming  effect  on  the  number  of  positrons  that  survive. 
Therefore  a  small  refinement  in  the  density  profiles  of  these  models  to  produce  a  better  fit 
to  the  light  curve  would  greatly  affect  the  results  of  my  model. 

This  strong  dependence  on  the  density  suggests  that  Type  I  supemovae  with  their 
larger  explosion  energies  (and  thus  lower  densities  at  a  given  time)  are  good  candidates 
for  the  source  of  positrons.  All  that  is  needed  is  a  nearby  Type  I  supernova  to  improve 
theoretical  models. 

In  my  computer  model  I  used  positron  injection  times  ranging  from  5  to  30  years. 
Figure  14  shows  that  a  peak  level  of  production  was  reached  within  that  time  period  for 
each  shell,  ie.  the  injection  time  that  produced  the  most  "survivors"  from  each  shell  was 
between  5  and  30  years.  Thus  I  feel  that  my  models  adequately  describe  the  production 
of  Galactic  positrons  from  Type  II  supernovae  given  their  somewhat  limiting 
assumptions.  More  complicated  hydrodynamic  models  which  take  into  account  the 
structure  of  the  magnetic  fields  and  the  details  of  the  radial  mixing  may  yield  more  precise 
results  but  they  should  not  negate  my  final  result. 


Figure  14:  This  figure  shows  the  injection  time  when  the  largest  contribution  was 
made  to  the  "survivors"  from  that  shell.  Note  that  the  shells  of  lowest 
density  made  their  largest  contributions  between  6  and  12  years. 
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Appendix  A:  Comparison  of  Analytic  and  Mom 


Program:  fdecay.c  Last  revision:  1/18/89 
Time  Monte  Carlo  Survival  Prob. 


9.460800e+08 

9.776160e+08 

1.009152e+09 

1.040688e+09 

1.072224e+09 

1.103760e+09 

1.135296e+09 

1.166832e+09 

1.198368e+09 

1.229904e+09 

1.261440e+09 

1.292976e+09 

1.3245 12e+09 

1.356048e+09 

1.387584e+09 

1.419 120e+09 

1.450656e+09 

1.482192e+09 

1.513728e+09 

1.545264e+09 


1.040000e-03 

1.600000e-03 

2.520000e-03 

3.880000e-03 

5.460000e-03 

7.460000e-03 

1.014000e-02 

1.302000e-02 

1.658000e-02 

2.120000e-02 

2.698000e-02 

3.310000e-02 

4.022000e-02 

4.758000e-02 

5.568000e-02 

6.472000e-02 

7.426000e-02 

8.482000e-02 

9.612000e-02 

1.080000e-01 


Carlo  Survival  Routines 


Analytic  Survival  Prob. 


1.035319e-03 

1. 658541  e-03 

2.543862e-03 

3.755 176e-03 

5.358753e-03 

7.420789e-03 

1.000504e-02 

1.317069e-02 

1.697060e-02 

2.144997e-02 

2.664539e-02 

3.258440e-02 

3.928528e-02 

4.675734e-02 

5.500124e-02 

6.400976e-02 

7.376845e-02 

8.425654e-02 

9.5447 82e-02 

1.0731 15e-01 


it  ri  in  ■  wiiii  or 


INTT.  COND. 


SHELL  -  1 

S  MAX  * 

SHELLS 

SHELL-1 

^  YES 

TIME -START 

S  STOP 

TIME 

TIME-STEP 

^  YES 

ENERGY-ST ART 

S  STOP 

ENERGY 

ENERGY-STEP 

ENERGY!  THERM 


YES 


CALCULATE  dE /dx.  dE.  dt 


FIND  TIME  OF  THERM. 


CALCULATE  MEAN  THERM 


TE  SURVIVAL  PROB 
mi  NUMBER  of  SURVIVORS 


CALCULATE  «  FROM  S 


CALCULATE  *  FROM  ALL 
SHELLS 
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/***•* . .  LiFETTMEW.C  4/19/89 

//include  <stdio.h> 

#include  <malh.h> 

#define  CUT_OFF_ENERGY  2.6e-4 
#defme  INniAL.DENSITY  1.08el0 
graph**/ 

#define  DsTTTAL_TIME  7.43e6 
seconds**/ 

#define  RAD_RECOM_RATE  2.0e-12 
#defme  COBALT_SOURCE  1.4987e54 
atoms**/ 

//define  DECAY_RATE  1.03785e-07 
//define  ION_POT  2.99e4 
#define  ZJFREE  1.5 
electrons/atom**/ 

#defme  Z_BOUND  24.5 

#defme  IN_FUGHT_FRACTION  0.0 

flight**/ 

//define  REST_MASS  0.511 
//define  ENDJTTME  3.473495e9 

FILE  ‘point; 

/**********************************************< 


/*******************************  SUBROUTINES  ***' 
/********************************************— *** 

/••••••Ionization  Losses  using  the  handbook  formula******/ 
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/**  inMeV**/ 

/**atoms/cc,  from  Woosley's 

/**86  days  expressed  in 

/**Lrom  Bussard's  paper**/ 

/** 0.07  solar  masses  converted  to 

/••1/111.52  days**/ 

/**11.5*Z  inMev**/ 

/••assuming  1 5  free 


/••Number  that  annihilate  in 


/••when  tau  is  100  years**/ 


7 

**/ 
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double  ionization/energy,  density,  betasqrd) 
double  energy,  density,  betasqrd; 

{ 

double  ratio,  parti,  pan2,  part3,  dEdxl ; 

ratio  *  pow(energy/REST_MASS,  0-5); 

parti  =  5.1e-25*(density*Z_BOUND)/betasqrd; 

part2  *  pow(betasqrd,  0.5)*(energy+REST_MASS)*ratio/ION_POT; 

part3  *  0.5*betasqrd; 

dEdxl  =  partl*(log(pan2)-part3); 

return  (dEdxl); 

} 


/*’ 
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/••••••Plasma  Losses  using  Jackson's  formula******/ 

double  plasma(  energy,  density,  betasqrd) 
double  energy,  density,  betasqrd; 

{ 

double  gamma,  part4,  part5,  part6,  dEdx2; 

gamma  *  pow((l -betasqrd),  -0.5); 

part4  *  5.1e-25*(density*Z_FREE)/betasqrd; 

part5  »  53e5*pow (betasqrd,  03)*pow( (density *Z_FREE),  -0.5); 

partfi  m  3.86e-ll*pow((2.0/(gamma-1.0)),  0.5); 

dEdx2  ■  part4*log(0.764*part5/part6); 

return  (dEdx2); 

} 

/******************************************************************************/ 


/******Racliatjve  Capture  of  Positrons******/ 

double  decay(meamherm,  initdensity) 
double  meantherm,  initdensity; 

{ 

double  fraction,  stop,  taunot,  timesqrd,  const; 

timesqrd  =  pow(mear therm,  -2.0); 
stop  =  pow(END_TTME,  -2.0); 

taunot  =  l/(initdensity*Z_FREE*  R  AD_RECOM_RATE*po  w(IN  11IAL_TTME,  3 .0)); 
const  =  pow((taunot*2),  -1.0); 
fraction  =  exp(const*(stop  -  timesqrd)); 

retum(fraction); 

) 


/* 
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/••••••Number  of  Positrons  Injected  at  Starttime******/ 

double  injection(starttime,  step,  shellsource) 
double  starttime,  step,  shellsource; 

( 

double  start,  finish,  num; 

num=shellsource*(exp(-(DECAY_RATE*(stantime-(step/2))))  - 
exp(-(DECAY_RATE*(starttime+<step/2))))); 

retum(num); 

} 


/******************************************************************************/ 


/*••••**••*******•*••******•****  MAIN  PROGRAM  ******••**•****••••*••******•*' 
/******************************************************************************/ 


mainO 

{ 

double  energy,  startEn,  time,  starttime,  step; 

double  density,  scale,  mean  therm; 

double  dEdx,  dEdxl,  dEdx2,  dE,  aveenergy; 

double  velocity,  dt,  total,  survivors,  numinjected,  percentleft; 

double  avebeta,  betasqrd; 

double  store[10][2],  distrib[10J.  shell[15][2]; 

double  initdensity,  shellsource,  grand  to  til,  sum; 

int  x,yA»; 

point*fopen("dataW","w"); 


fprintf (po  in  t,’VNn"); 

fprintf  (point,"  PROGRAM:  LifetimeW.c  LAST  REVISION:  4/19/89Wn"); 

fprintf (point,"  Start:  S  yean  Finish:  30  yean  Step-  1  yearin'm'n*); 

fprintf  (point," - Wn"); 

distrib[0]»  0.061; 
disoibtl]-  0.117; 
distrib[2]>  0.1S2; 
distrib(3]-  0.165; 
distrib[4]>  0.158; 


distrib[5]=  0.142; 
distrib[6)=  0.104; 
distrib{7]=  0.066; 
discrib(8]=  0.030; 
distrib[9]=  0.005; 

shell[0][0]=  1.08el0; 
shell[l][0)=  5.38el0; 
shell[2][0]=  3.76ell; 
shell[3][0]=  5.38el0; 
sheU(4][0]=  3.76el0; 
sheU[5][0]=  2.37el0; 
sheU[6][0]=  2.37el0; 
sheU(7]{0]=  2.16eI0; 
sheU[8][0]=  1.61el0; 
sheU[9][0]=  3.23e9; 
sheU[10][0]=2.16e9; 
sheU[ll][0]=2.69e9; 
sheU[12][0]=3.23e9; 
shell[13][0]=  3.23e9; 
shell[14][0]=  3.23e9; 

shell[0]l  1]=.  2.246e53; 
shell[l](l]=  2.332e53; 
shcll[2][l]=  1.925 e53; 
sheU[3][l]=  1.562e53; 
shell[4][l]=  1.305e53; 
shell[5][l]=  1.070e53; 
shell[6][l]=  8.770e52; 
shell[7][l]«  7.273 e5  2; 
sheU[8][l]=  5.773 e5  2; 
shell[9][l]=  4.706e52; 
shell[10][l]=  4.064e52 
shell[ll][l]=  3.209e52 
shell[12][l]=  2.567e52 
shell[13][l]=  1.925e52 
shell(14J[l]=  1.7llc52 


grandtotal »  0; 
step  =  3.1536e7; 
z  =  0; 

r*”***””*”*”*”**”*»»START  OF  SHELL  LOOP*********************.*****/ 

for  (z=0;  z<  14.5;  z*z+1.0) 

( 

total  *0; 

initdensity  *  shell[z][0]; 
shellsource  *  shell[z][l]/3; 

fprintf(point,“  Shell  number.  %2.0f  Initial  density:  %53e  Amount  of  cobalt:  %53^n",  z+1.0, 
initdensity,  shellsource); 

fprintfCpoint,"  INJECTION  TIME  TIME  OF  THERM.  NUMBER  OF  POSITRONS  PROB.OF 
SURVIVALS); 

fprintf(point," - - - — - - . - -\nH); 


r 


START  OF  TIME  LOOP*' 
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for  (starttime=  1.5768e8;  slant ime<=  9.47e8;  starttime=starttime  +  step) 

( 

x=0; 

y=0; 

meantherm=0; 

/•• . •••••START  OF  ENERGY  LOOP*************** 

for  (startEn=  0.073;  surtEn<=  1.46;  stanEn=startEn  +  0.146) 

( 

time=s  tan  time; 
energy=startEn; 

LOOP*******************' 

while  (energy  >  CUT_OFF_ENERGY ) 

( 

scale=(5.0el0/3.92e21)*pow(time,3.0)*shell[0][0]/sheU[z][0]; 
density=initdensity*pow(INlTlAL_TTME,3.0)/pow(time,3 .0); 

betasqrd=l-pow(((energy/REST_MASS)+1.0),  -2.0); 
dEdxl=ionization(energy,  density,  betasqrd) 
dEdx2=plasma(energy,  density,  betasqrd); 
dEdx=dEdxl+dEdx2; 
dE=scale*dEdx; 
if(energy-dE>0) 

{ 

aveencrgy=energy-(dE/2); 

avebeta=po w(  1  -po w(((aveenergy/REST_M ASS)+1 .0),  -2.0),  0.5); 

velocity=3el0*avebeta; 

dt=scaIe/VeIocicy; 

/••••••••••InCTement  Time  and  Energy**********/ 

energy=energy-dE; 
times  time+dt; 

) 

else  goto  end; 

} 


/.*.*.«****.***«*««**«*. END  of  WHILE  LOOP*** 
/••••••♦••♦Calculate  and  Store  Lifetime**********/ 

end:; 

store[x-H-][0]*time; 

) 

/.........•...............END  OF  ENERGY  LOOP*' 

/.*.**•.*. ‘Calculate  Mean  Lifetime**********/ 

while  (y  <  x) 

{ 

«ore[y][l]«stoie[y][0]*diitrib[y]; 
mean  therm»mean  therm -«tore(y]  ( 1  ] ; 

y*y+i: 

} 


r 


'Calculate  How  Many  Survived' 
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percentleft  =  decay(meanthetm,  initdensity); 
numinjected  =  injection(staraime,  step,  shellsouree); 

survivors  =  percentleft*(numinjected*(l-IN_FLIGHT_FRACnON)); 
total  =  total  +  survivors; 

fprintf(point,"  %t  %c  %e\n",starttimejneamhermjiuminjected,percentleft); 

} 

****•••***•*•***••*•**•*••  OF  TIME  LOOP****************************/ 

fprintf(point,'\iNn'');  fprintf(point,"Number  surviving  from  shell  %1.0f:  %eNn'>n",  z+1.0,  total); 
fprintf(point," - \nV); 

grandtotal  =  grandtotal  +  total; 


/***********«**********«*****£NQ  OF  SHELL  LOOP****************************/ 

sum=0; 

for(s=0;s<14.5;s=s+l) 

l 

sum=sum+<shell[sn  1  ]/5); 

) 

fprintf(point,"The  total  number  surviving  =  %e  out  of  %e.\nW,grandtotal,sum); 
fprintf(point,'The  survival  fraction  is  %e.",  (grandtotal/sum)); 

fclose(point); 


1 

/****'****************•****“**********“***************“*********************/ 

/*********»****************«****TJ{£  £J^O**************************************** 

. . . / 
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